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We compare recent numerical results, obtained within a "helical Killing vector" (HKV) approach, 
on circular orbits of corotating binary black holes to the analytical predictions made by the effective 
one body (EOB) method (which has been recently extended to the case of spinning bodies). On 
the scale of the differences between the results obtained by different numerical methods, we find 
good agreement between numerical data and analytical predictions for several invariant functions 
describing the dynamical properties of circular orbits. This agreement is robust against the post- 
Newtonian accuracy used for the analytical estimates, as well as under choices of resummation 
method for the EOB "effective potential", and gets better as one uses a higher post-Newtonian 
accuracy. These findings open the way to a significant "merging" of analytical and numerical 
methods, i.e. to matching an EOB-based analytical description of the (early and late) inspiral, up 
to the beginning of the plunge, to a numerical description of the plunge and merger. We illustrate 
also the "flexibility" of the EOB approach, i.e. the possibility of determining some "best fit" values 
for the analytical parameters by comparison with numerical data. 

PACS numbers: 04.30.Db, 04.25. Nx, 04.25. Dm, 04.70. Bw, 97.60.Lf, 97.80.-d 



I. INTRODUCTION 



Binary black holes are the most promising candidate sources for the LIGO / VIRGO / GEO600 / TAMA . . . 
network of ground based gravitational wave (GW) interferometric detectors jl], |2|, ||, p|. Signal to noise ratio 
estimates suggest that the first detections will concern black hole binaries of total mass £ 25Af Q , and show that 
the most "useful" part of the gravitational waveform is emitted in the last ~ 5 orbits of the inspiral, and during the 
plunge taking place after crossing the last stable (circular) orbit (LSO) S. This makes it urgent to have reliable 
methods allowing one to model the last orbits of binary black holes. 

Recently, Ref. Q has suggested a new method to tackle the dynamics and GW emission from the last orbits of 
binary black holes. The basic idea of |(| was to extend, by using suitable resummation methods, the validity of the 
pcrturbative ( "post-Newtonian" ) analytical calculations of the equations of motion 0, ||, ||, and GW radiation 
pH [l2[ so as to be able to describe not only the last orbits before the LSO, but also the transition to the plunge. 
The main philosophy of || was to push analytical methods to their limits so as to use numerical methods only to 
describe the plunge, which was estimated (when the masses are comparable) to last about 60% of one orbit. To 
implement this philosophy, one needs to be able to construct numerical initial data that match the dynamical initial 
data given by the method of [|| at the start of the plunge. This is a non trivial task because, up to very recently, there 
was a significant discrepancy between analytical || |l3|, |l4|, [l5j and numerical |H| [D], |l^| estimates of the dynamical 
characteristics of the orbits near the LSO. For instance, the binding energy at the LSO, eLso — ^lso/(™i + m 2) — L 
is analytically estimated to be around ef^o — —1-67% [Oj (third post-Newtonian effective-one-body estimate with 
uj s = ||), while two different numerical estimates |l6|, |ll| gave e£g™ — —2.3%, i.e. a significantly more bound (by 
38%) LSO. The discrepancy is even more striking if one considers the LSO orbital period: the analytical estimate is 
T£g — 71.2 (mi + 7712) [Q, which is twice longer than the numerical one T^g™ — 35 (mi + 7712) [ [f8| - [By contrast, 
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the dispersion among analytical estimates || [l3|, [l4|, |T^] , as well as the dispersion among numerical ones |l6|, [I?], [T^] 
is significantly smaller than the difference between analytical and numerical results. See |19| for a review of analytical 
methods, and see below for a brief discussion of our choice of analytical estimates.] This discrepancy diminishes the 
physical relevance of a recent attempt [^o) at fulfilling the proposal of Ref. |(|, i.e. to start a full numerical calculation 
of the plunge just after the crossing of the LSO. 

Very recentlv, a new numerical approach to the circular orbits of binary black holes has been set up [plf and 
implemented [£2| (see also Ref. for a discussion and an extension to any black hole rotation state). Contrary to 
the previous numerical approaches jl^, [l8| , which were formulated in terms of the initial value problem of general 
relativity and dealt therefore only with the four constraint equations on a 3-dimensional spacelike hypersurface, the 
new approach |2^] deals with a full 4-dimensional spacetime (see Ref. [^3| for the link with York's conformal thin- 
sandwich formalism |^4|, ^5|). The basic assumption is that this spacetime is endowed with a helical Killing vector 
pq| , which amounts to neglecting the gravitational radiation reaction, i.e. to considering that the orbits are exactly 
circular. In addition, Refs. |22| ] assume, as a simplifying approximation, that the spatial 3-metric is conformally 
flat (see Sec. IV. C of Ref. [p7| for a discussion). The slicing is also chosen to be maximal (i.e. K = 0). The spacetime 
manifold is chosen to have the topology of the real line times the two-sheeted Misner-Lindquist manifold. The two 
sheets are assumed to be isometric (with respect to the full 4-metric). Under the above assumptions, the problem 
amounts to solving five of the ten Einstein equations. This is one more that the previous numerical approaches 
0, 0, E8l Hereafter we call the new method (2lJ Q the helical Killing vector (HKV) approach, and previous 
methods |16|, |l7| , |Tg|| IVP ones (for Initial Value Problem). Note that the HKV approach has also been employed for 
binary neutron stars (see e.g. p8[ and references therein). 

The first numerical results obtained by the HKV method [^2| for the characteristics of the LSO indicate a much 
better agreement with analytical estimates than those obtained with the IVP method |l7|, [l8|]. However, the 
comparison, done in p2| , used analytical results valid only for non spinning black holes, while the numerical approach 
of [ pT| concerns corotating binary black holes. Moreover, the comparison made in p2| was limited to the sole LSO, 
without looking at the behavior of the other circular orbits. It is known that the presence of even small additional 
interactions (especially repulsive ones) can have an important effect on the dynamical characteristics of the orbits 
around the LSO. It is therefore a priori necessary to include spin-orbit effects (which are repulsive in the corotating 
case) when doing the analytical/numerical comparison. Fortunately, a recent analytical work [^9| has shown how to 
generalize the effective one body (EOB) approach of Q by including, to lowest order, spin-orbit and spin-spin effects. 
The purpose of the present work is to take into account in detail the effect of the spin interactions in corotating binary 
black holes. Another important feature of this work is that we shall compare the predictions of the EOB method 
H ^9) and the results of the new numerical approach of |2l], , not only for the characteristics of the LSO (to 
which most authors limit their considerations), but also for all other circular orbits. We shall do that by comparing 
several physically important (and coordinate-invariant) functions, such as the binding energy as a function of angular 
velocity. Note that, in this work, we loosely call "LSO", for a corotating system, the maximal binding configuration 
along a corotating sequence. For HKV sequences, it has been shown that this minimal energy configuration marks 
the onset of an orbital instability j2?j . The fact that such corotating sequences are (probably) not realized in actual 
coalescing systems 1 , and that the corotating "LSO" signals only a secular instability, is not important for the present 
work whose main aim is to show how one can relate the numerical and analytical descriptions of the same configuration 
involving spinning black holes. 

This paper is organized as follows. Section II discusses the arguments selecting the EOB method as the current best 
analytical method for tackling the last orbits of binary black holes and summarizes the EOB approach to spinning 
black holes. Section III applies the EOB formalism to corotating black holes on circular orbits and shows how to 
compare the analytical results to the recent numerical computation of a sequence of corotating binary black holes. 
Section IV discusses the meaning of our results. 

II. EFFECTIVE ONE-BODY APPROACH TO SPINNING BINARY BLACK HOLES 

A. Formulation 

Before summarizing the effective one-body (EOB) formalism, let us recall why we consider it as the best current 
analytical formalism for tackling the dynamics of the last orbits of binary black holes. The first point is that the 



Note, however, the claim of [ pp| that the effective viscosity of black holes might be large enough to trigger a significant angular-momentum 
transfer, potentially able to ensure tidal locking at small separations. 



3 



second work in Ref. || has shown that it is crucial, if one wishes to keep a good overlap with the expected real 
signals, to generate GW templates which go beyond the "adiabatic approximation" and which use accurate equations 
of motion, including radiation reaction, to describe the smooth transition (taking place around the LSO) between 
inspiral motion and plunge. However, the most accurate equations of motion [|, are a priori given as very 
complicated expressions, containing hundreds of terms. These equations of motion are essentially power series in a 
formal "small parameter" e ~ GM/c 2 r ~ v 2 /c 2 ["post- Newtonian" (PN) expansion]. The problem is that the "small" 
PN expansion parameter e is not numerically very small near the LSO. It was found quite early that the use of 
straightforward PN-expanded equations of motion (as in is unsatisfactory because the PN-expanded equations of 
motion converge so slowly that, at low orders,they can lead to a behavior which is qualitatively different from what 
one expects. In particular, it was shown in |32] that the PN-expanded (harmonic-coordinates) equations of motion 
admit an LSO at the second post-Newtonian (2PN) level, but admit no LSO at both the first post-Newtonian (1PN) 
and the third post-Newtonian (3PN) levels. [ "LSO" is taken here in the sense of a critical point in the dynamical 
analysis of circular orbits.] By contrast, a study of the LSO in terms of the PN-expanded Hamiltonian |33j has found 
that the PN-expanded Hamiltonian (in various coordinates) admits no LSO at the 2PN level, but admit LSO's at the 
1PN and 3PN levels. This unstable behavior is a sign of the unreliability of non-resummed PN approaches. 

As the recent work [HJ has succeeded in deriving complete 3PN equations of motion, it would be unacceptable 
to "spoil" the precious information they contain by using them in such an inappropriate way. This led Ref. |p2fl to 
propose to improve the situation by using a partial resummation approach. Specifically, they introduced an "hybrid" 
approach in which one resums exactly the "Schwarzschild" terms in the (relative) equations of motion, while continuing 
to treat the ^-dependent terms (where v = m\m<il{m\ + 7712) 2 ) as straightforward PN expansions. This "hybrid" 
approach did improve the situation to some extent (at least at the 2PN accuracy). However, later work showed 
that the hybrid approach was neither robust, nor consistent. Ref. p3| showed that the predictions of the hybrid 
approach were robust neither under using Hamiltonian equations of motion instead of Newtonian-like ones, nor under 
a change of coordinate system. Ref. p| pointed out the inconsistency of the hybrid method by showing that the 
formal (un-resummed) "^-corrections" represent, in several cases, a very large (larger than 100%) modification of the 
corresponding resummed ^-independent terms. This led to a quest for new resummation approaches providing full 
dynamical equations of motion for binary systems, and being robust and consistent. The only such method we are 
aware of is the EOB approach jfj). We do not consider here methods such as the E- or e-methods of 13 and the 



j-method of Q] which were devised to investigate the location of the LSO. These methods cannot be used beyond 
the "adiabatic approximation" because they do not provide full dynamical equations of motion able to model the 
transition between inspiral and plunge. Let us note in this respect that the work of Q has shown that, in the case of 
comparable masses, the transition between inspiral and plunge is a gradual process in which the precise location of the 
"adiabatic" LSO is less important than a precise evaluation of radiation damping effects. Therefore it is finally not 
very important to focus on the "exact" location of the LSO. What is really important is to have a good description of 
the dynamics of the last ten orbits before the plunge. Moreover, even if we were to consider only the determination 
of the LSO, the comparisons made in jl4j have shown that the EOB method is the most robust in that it is the least 
sensitive to changes in the 3PN coefficients (see Fig. 2 in [Q). We also note in passing that |l3| had shown that 
the e-method was more robust (both in its dependence on v and under changes in higher PN coefficients), and more 
convergent, than the straightforward non resummed .E-method (recently used at 3PN in Jlj|); see Table I, Table VII 
and Fig. 5 in [||. The robustness of the EOB method is discussed in Refs. §, 0, §9). We shall give below further 
examples of the robustness of the EOB approach. 

The basic idea of the EOB method is to map the (complicated) relative dynamics (in the center of mass) of a 
two-body system onto the (drastically simpler) dynamics of an "effective" one-body system. Most of the (physically 
irrelevant) gauge-related complications of the two-body equations of motion get absorbed into the mapping between 
the two problems. The mapping is defined so as to preserve the adiabatic invariants (which, at the quantum level, are 
quantized in units of K). The energy mapping E c g — f(E rea x) between the real two-body (center of mass) energy E lea i 
and the effective one-body energy E c g is determined by the matching between the two problems. It is remarkably 
found that the energy map / is extremely simple, and independent of the PN accuracy considered (we set c = 1): 

E eS _ £wi -m\- m\ 



fi 2mi m2 



(1) 



Here, mi and ni2 denote the masses of the two bodies, while /i denotes the reduced mass fi — mims/M, with 
M = mi + 1^2 denoting the total mass. Let us note also the definition of the symmetric mass ratio v = [ijM = 
mi m2/(mi + m-i) 1 (0 < v < i). The inverse of the energy map (^) is 




E ieal = MJl + 2u(^-^). (2) 
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The effective energy E e g is given by the effective Hamiltonian H e g(x,p, Si, Sy. Here, x and p are the positions 
and momenta of the effective one-body (describing the relative motion) and Si and S2 are the two independent 
spin vectors of the two bodies. Ref. |(| considered only the 2PN approximation and the spin-independent case. The 
3PN, spin- independent effective Hamiltonian was determined in [Q. Ref. |^9| recently showed how to add the spin 
interactions in the EOB framework. The most complete version of the 3PN spin-dependent effective Hamiltonian 
is given by Eq. (2.56) of [^9|. For simplicity, and because spin-spin effects will turn out to be quite small, we shall 
consider the less complete, but simpler, effective Hamiltonian defined by Eq. (2.26) of |^9| with the effective metric 
defined by the "deformed Kerr metric" of Section IIC there. This Hamiltonian takes into account (at the lowest 
PN approximation) the full spin-orbit interaction, and incorporates most of the spin-spin interaction. [In the case of 
parallel spins the ratio between the spin-spin interaction included in the simplified H c g and the more complete one 
H' cS is (|) .] The Kerr parameter of the effective Kerr metric describing spin-orbit and spin-spin interactions is given 
by 

Ma = S e s = <J\ Si + a 2 S 2 , (3) 

where 01 = 1 + 3m2/(4mi), 02 = 1 + 3mi/(4m2). 

The EOB approach can deal with the most general case where the spin vectors Si, S2 are arbitrarily oriented. 
In fact, one of the great advantages of the EOB approach is that, being Hamiltonian-based, it provides a full set 
of evolution equations for all the dynamical variables: x, p, Si and S2. See Eqs. (3.1)-(3.4) of ]2£|. Here, we are 
interested in comparing the EOB predictions with the numerical results of p2f which are restricted to the simple 
case where Si and S2 are aligned with the orbital angular momentum L. In this case, it was shown in |29[ ] that, 
for a given effective energy E c g, and a given angular momentum L, the (effective) radius r = \x\ of the circular 
orbits is obtained by solving the two equations, R(r, E e Q, L) = R(r, E e g, L) = 0, where R is a certain effective 
radial potential. [We use the fact that the Carter-like constant Q vanishes for the "equatorial" orbits that we are 
considering.] Actually, it is more convenient to work with the inverse radial variable u = M/r (we set G = 1) and 
with the corresponding w-potential U{u) — yT 2 r -4 R(r). Let us also introduce the following dimensionless variables 
(note that L = L z for the orbits we consider and that, by definition, L = L e a — L re ai) 

s _ E roa i - L _ a 

E — , L = , a = — , (4) 

H ' nM M w 

and 

E=^ff, L^ L ~ aEcS ^L-aE. (5) 

The quantity a entering these equations is the modulus of the effective Kerr parameter defined by Eq. (^). 
In terms of these quantities the (scaled) u-potential reads 

U(u) = (E-aLu 2 ) 2 -A{u)(l + L 2 u 2 ), (6) 

where the function A(u) = A(u) + a 2 u 2 will be discussed below. 

The sequence of circular orbits is defined by the solutions of the two equations 

U(u,E,L) = 0=^U(u,E,L). (7) 

The effective Hamiltonian E = H e s(u,L) would be obtained by solving U(u,E,L — aE) = which is a quadratic 
equation in E. However, the intrinsic simplicity of the EOB approach is more visible if we use the variable L = L — a E 
instead of L = L/fxM. Indeed, the solution of U(u, E, L) = reads 



E = aLu 2 + \/ ' A(u){l + L 2 u 2 ) . (8) 

This result can be viewed as a simple modification by spin-orbit (a L u 2 term) and spin-spin effects (a 2 u 2 contribution 
to A(u) = A(u) +a 2 u 2 ) of the radial potential for circular orbits in a spherically symmetric gravitational potential 



(such as a Schwarzschild metric, or the effective metric in absence of spins): Eq(u, L) = y A{u)(\ + L 2 u 2 ). In such a 
situation A(u) would simply be (minus) the time-time component of the (effective) metric (1 — 2u in the Schwarzschild 
case). A nice feature of the EOB approach is that, in the spin-independent case, the dynamics of circular orbits is fully 
encoded in only one function: the function A(u) = — g^Q (u), with u = M/r. [As discussed in M, [l4|, |29|, the dynamics 
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of non circular orbits involves two more functions: D(u) and Q 4 (u,p).] It is truly remarkable that the hundreds of 
complicated contributions entering the PN-expanded two-body equations of motion get drastically compactificd in so 
few functions. In particular, for circular orbits, the essential physical information of the PN-expanded dynamics is 
contained in the effective metric coefficient 

A{u) = 1 - 2u + 2v u 3 + a 4 (u) u 4 + 0(u 5 ) . (9) 

Note that there are no 1PN (oc u 2 ) contribution in A(u), and that the 2PN effects amount to the very simple term 
§). The 3PN contribution oc u 4 was found jl4| to have the coefficient 

/94 41 \ 
a i (v)={—-—ir 2 + 2uj s \ v. (10) 

Again, remarkable simplifications happened in the calculation of the 3PN effective metric. The 3PN coefficient a^{v) 
could, a priori, have involved contributions proportional to v 2 and v 3 . All such contributions canceled to leave a 
simple final result 0.4 (V) oc v. The parameter uj s entering Eq. (|l0|) ( which had been left ambiguous by the first 3PN 
calculations) has been recently determined to be simply oj s = 0. This yields a numerical value for the general 
rclativistic prediction for the 3PN coefficient equal to: 

a<f R (z/) ~ I8.6881/ - 4.672 (4i/) . (11) 



B. Representation of the basic EOB function A(u) 

The EOB method a priori determines (from the original, fully PN-expanded dynamics) only the PN-expansion of 
the basic function A(u), or of the combination 

A(u) = A(u)+a 2 u 2 = 1 - 2u + a 2 u 2 + 2vu 3 + a 4 (v) u 4 + 0{u 5 ) , (12) 

entering the effective radial potential (0). 

In this respect let us emphasize the various senses in which the EOB approach is a "resummation" of the original, 
fully PN-expanded equations of motion. A first, crucial "resummation" feature of the EOB approach is the mapping 
of the original dynamical variables q 1 , p x , q 2 , P2 onto the effective variables x, p. This mapping 0, [l4j is given 
as a complicated PN expansion, and provides a first level of simplification and compactification of the dynamics. A 
second feature consists in the fact that the real Hamiltonian is given by an iterated square-root. For instance, in the 
spin-independent case (p = pj \i) 



H lea ,i(x,p) = M A 




( Mn l)(n-p) 2 + Q i {p)\ I. (13) 



This iterated square-root structure is, by itself, a partial resummation because the fully PN-expanded Hamiltonian, 
which has the form i? roa i (p) = M [1 + v ( | p 2 + ^pj + • • • + c 8 p 8 + ■ ■ ■ + O (p 10 )] , would be obtained by re-expanding 

all this square-root structure in powers of p 2 and M/r. Finally, on top of these two partial resummations, it is very 
natural, within the EOB approach, to further resum -ff rea i by replacing the PN-expanded results for A(u), Eq. (0), 
and/or A(u), Eq. ([l2]), by some better behaved expressions. Indeed, as we recalled above, in the spin-independent 
case A(u) plays the basic role of replacing — 3oo hwarz = 1 ~ 2M/r in the radial potential determining the circular 
orbits. 

It is then natural to require that the "exact" effective A(u) should have the same qualitative behavior as 1 — 2u, 
i.e. a simple zero at a finite value of u. More generally, in the spin-dependent case A(u) is a ^-deformation of the 
Kerr function r~ 2 A 4 (r) = r~ 2 (r 2 — 2 Mr + a 2 ) which determines the location of the horizon. It is therefore natural to 
require that the effective A(u) have, like 1 — 2u + a 2 u 2 , a simple zero at a finite value of u. To ensure such a qualitative 
behavior it is natural to replace the straightforward Taylor expansion giving A(u) by a suitable Pade approximant 
(chosen so as to robustly imply the presence of a simple zero) , namely |29J 

A(u,a 2 ) = P3 1 [1 - 2u + a 2 u 2 + 2v u 3 + a 4 (v) u 4 ] , (14) 

where denotes a Pade of the N n / D m type (where the indices denote the degrees of the numerator and denominator). 
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Various representations of A(u) 
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FIG. 1: Various representations of the function A(u,a 2 ). The right figure is a enlargement of the LSO region. 
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FIG. 2: Various Pade approximants of the function A(u,a 2 ). The right figure is a enlargement of the LSO region. 



As a test of the robustness of the EOB predictions we shall also consider a more brutal way of ensuring that A(u) 
has a simple zero at a finite value of u: it consists in denning 



A (u,a 2 ) = (1 - 2u)(l + 2vu z + (a 4 (» +4^)u 4 ) +a 2 it 2 



(15) 



Fig. |l| compares the various representations of the function A(u) : the straightforward Taylor scries approximant jl^ ) 
(truncated after the 3PN term oc it 4 ), the Pade approximant ( |l4| ) and the "factorized Taylor" approximant (|l5|). They 
are represented both in the limit a — 0, and for the value a = 0.15 which roughly corresponds to the value of a, at 
the LSO, in the effective one-body formalism discussed in the next Section. Note that the various representations 
start differing significantly from each other only for u ~ 0.4, i.e. for values about twice larger than the value of u 
corresponding to the LSO, i.e. ulso — 0.22. However, this does not mean that the resuming of A(u) and/or the 
addition of the spin-dependent interactions has no significant effect on the characteristics of the LSO. Indeed, as 
the usual (non corotating) LSO corresponds to an inflection point in the effective radial potential H c f[ (r, i, S a ), any 
modification of the radial dependence of the Hamiltonian (e.g. by a change in A(u) or the addition of the linear 
spin-orbit term aLu 2 in Eq. (||)) has a strong effect on the location of the LSO. 

On the other hand, Fig . ^| compares different variants of the Pade approximants that one could use. In the text, 
we have followed Ref. |2!| in focusing on the specific definition ( |l4| ) of the function A(u). However, one could have 
considered 



A'(u,a 2 ) = P£[l - 2u + 2vu A + a 4 (is)u 4 } +a 2 u 2 



(16) 



Another possibility arises when one uses the great "flexibility" of the EOB approach. As discussed in || g9| one can 
think of the EOB formalism as a multi-parameter analytical formalism, where only a fraction of the parameters is 
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currently known. For instance, there are further terms in the expansions (||) and (|l2|), say + a§(i>) u 5 + a${v) u 6 + ■ ■ ■ . 
The parameters a^lv), o,q{v), . . . are not known at present. However, they can be introduced as fitting parameters and 
adjusted so as to reproduce other information one has about the exact results. For instance, if we had reason to trust 
some numerical results, one could search for the optimal values of 04(1/), a§(v), . . . that best fit the numerical data. 
We shall give an example of this below. A minimal requirement is to test the "robustness" of the EOB formalism 
under adding some reasonable values for 125(2/), .... See Ref. pB| for a study of this robustness. Here, we just wish 
to illustrate a minimal aspect of this robustness: when adding the next term a§(y) u 5 to A(u), the natural Pade 
approximant to use becomes 

A'"{u, a 2 ) = Pl[l - 2u + a 2 u 2 + 2vu z + a 4 (is) u 4 + a b (v) u 5 } . (17) 

If we take the 05 — > limit of the Pade approximant (|l^) we do not recover (14). Indeed, the P3 1 approximant 
( |l4| ) is uniquely defined by the requirement of being oc N1/D3 and of matching (12) up to a^u 4 . This means that 
F3 effectively makes an "educated guess" about an infinite number of additional contributions as u 5 + a§ u 6 + • ■ • , 
which are plausible continuations of the known terms 1 — 2u + ■ ■ ■ + 04 u 4 . In particular, the Taylor expansion of P3 1 
must contain a non-zero contribution 05 u 5 . This explains why it does not agree with the 05 = limit of (^|). It is 
therefore interesting to compare ( |l7| ) to ( |l4| ) to measure the robustness of Pade approximants against knowledge or 
lack of knowledge of 05. In fact, Fig. || shows that the various ways we have just described of defining some Pade 
approximants of A(u) give extremely close results. Much closer than the variants represented in Fig. [j]. This is an 
illustration of the generic robustness of Pade resummation. 



III. COMPARING ANALYTICAL AND NUMERICAL RESULTS FOR COROTATING BINARY BLACK 

HOLES 

A. Corotating configurations in the EOB framework 

Let us now tackle the central problems of this paper: (i) to compute the coordinate-invariant functions characterizing 
the dynamics of adiabatic sequences of corotating binary black holes on circular orbits, and (ii) to compare the (EOB) 
analytical predictions with the numerical results of Refs. pl| , p2[ . The main subtlety in this study comes from the 
corotation assumption. This assumption means that the spins of the black holes vary along the sequence, and in 
fact increase as the radial distance r decreases. However, any point in the sequence, i.e. any circular orbit, must be 
determined (in the EOB approach) by imposing that 

„ . dH rcs x{r,p r ,L,Si,S 2 ) . . 

dr ' ( 8) 

„ . dH iesd (r,p r ,L,Si,S 2 ) 

0=Pr = ^ , (19) 

dp r 

and p r = 0, where the derivatives are all considered for fixed values of L, Si and S2. Eq. ([l9]) is automatically satisfied 
by p r — 0, and Eq. ( |l8|) (with p r = 0) then leads to Eqs. (0) above. Eq. ( |l8| ) determines r as a function of L, S\ and 
S2, or L as a function of r, Si and S2. It then remains to determine how Si and S2 vary as functions of L and r. Let 
us also note beforehand that, in view of Eq. (|l8|), the usual (irrotational) LSO corresponds to an inflection point (in r) 
of iJ roa i(r, L, S a ) considered for fixed values of L and S a (corresponding to the LSO), and to a minimum of H rea \ (i.e. 
a maximum binding configuration) along a sequence of fixed values of S a . In the corotating case, we define the LSO 
as the maximum binding configuration along a corotating sequence. The variation of S a with r along a corotating 
sequence changes the location of the minimum of iJ rca i with respect to the irrotational case. The variation of S a with 
r prevents this minimum to correspond to an inflection point of the function H tca \{r,L,S a ) considered for the fixed 
LSO values of L and S a . 

The variation of the magnitudes of the spin vectors along a sequence obliges us to complete the formalism of pjj 
by giving a prescription for the variation of the black hole masses mi or m 2 with S 2 or S 2 ,. [This issue did not show 
up in [p9[ because the EOB formalism guaranteed the conservation of S 2 and S 2 ,.] Though, in principle, there are 
subtleties linked to the exact physical meaning of the EOB spin vectors, at the level of approximation of Ref. |2j| 
it is clear that we can simply require that each gravitational (and inertial) mass m a (with a = 1,2) depends on S a 
essentially as in the Christodoulou mass formula Q 

^K"< 8 a ] = J (mr) 2 + 4^55 • (20) 
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Here, the irreducible mass m 1 ™ is an adiabatic invariant (constant along a corotating sequence), which should be 
related to the area of the horizon via A a = 167r(mJ t rr ) 2 . 

The simplest definition of a corotating sequence is that the spin angular velocities of the black holes fl a stay equal 
to the orbital angular velocity Qq: 

Oi = Q 2 = O . (21) 

The universal relations of (Hamiltonian) dynamics (such as 17 = dif/dt = + dH rca _\/ dp v with p v = L z = L) tell us 
that the angular velocities are obtained (for circular orbits with aligned spins) as 

dH iea x(r,L,Si,S 2 ) dH Ieal (r ) L,S 1 ,S 2 ) . . 

n ° = di ' n ° = as a ' (22) 

where H IC& i(r, L,S a ) is the full (real) Hamiltonian, 

-ffrcai = H Ica ,i(m'", r,p r ,p v , S a ) , (23) 

considered for fixed values of the irreducible masses m 1 ™ , and for p r = (and p v = L z = L). Note that this means, in 
particular, that in expressions such as Eq. (^|) the masses mi, m 2 , as well as M = mi+m 2 and v = mi 1122/ (mi+m2) 2 , 
must be considered as functions of m 1 ™ and S a ■ The dependence of the masses on S a is crucial in the computation of 
the spin angular velocities Q a = dH Tea i/dS a . In fact, one can see that f2 a is the sum of two contributions: a "bare" 
contribution f^ arc = dm a /dS a which is the angular velocity of an isolated Kerr hole, and a contribution linked to 
the presence of the other hole, which represents a frame-dragging effect, which is automatically included in the EOB 
approach. [Note the consistency of the EOB approach under the replacement (pp|): this yields a new (numerically 
dominant) contribution to £! a = d ^g al which drops out of the spin-evolution equation (3.5) of Ref. |^9|. Note also 
that one would need to introduce extra (dissipative) terms in the spin evolution equation to enforce a realization of 
a corotating sequence in the EOB formalism.] 

Along any continuous sequence of circular orbits (not necessarily a corotating one) one has dH Tea \/dr = = 
dH Tea x/dp r (see Eqs. (|l8|), (|l9|)), so that we can write the variation of the energy as 

dE rca i = n dL + fti dSi + Q 2 dS 2 ■ (24) 

Note that, in the special case of a corotating sequence we have 

eLKrcai = n d J (25) 

where J = L + Si + S 2 is the conserved total angular momentum of the EOB approach (see Eq. (3.7) of ^9|), and 
where f2 is the common value of the angular velocities ( |2~]] ) . 

Let us, for a moment, shift from the analytical side to the numerical side. One of the nice advantages of the 
numerical approach of |22| is that it determines not only the global (Arnowitt-Deser-Misner; ADM) conserved 
quantities Madm and Jadm, but also the common angular velocity fix (which appears in the asymptotic boundary 
condition for the helical Killing vector i^d^). A corotating sequence is then numerically defined by imposing the 
satisfaction of gJ-Madm = &k (^Jadm- In view of Eq. p5]) , it is then fully consistent to identify the numerical and 
analytical quantities according to Madm = E Ica \, Q,k = O, Jadm = J- Ref. Q has also defined a total (numerical) 
irreducible mass M- lrr as the sum ^Ai/Wir + yM^>/167r, and has checked that it remains constant along a corotating 
sequence. It is therefore fully consistent to identify M- m with the analytical quantity m 1 1 rr + m 1 2 rr . [Once the constancy 
of Mj rr is numerically verified, the consideration of Mj rr for infinitely separated configurations, for which the areas A a 
unambiguously tend to their Kerr values in both formalisms, forces the identification M- m = m" r + m" 1 '.] 

The numerical results of [pjl give access to several invariant functions, such as Maj^m/M^ as a function of £Ik M; rr , 
and Jadm/M^j. as a function of VIkM^. However, let us mention that this normalization with respect to M- lrl is 
different from the one used in , where the global quantities are chosen so that Madm = 1 at the LSO (normalization 



with Mo in the language of 22 ). Now that we have shown that these quantities have unambiguous correspondents 
in the (EOB) analytical framework, it remains to complete the analytical determination of the corotating sequence. 
To do this we need to write explicitly the equations (0) and (21). Let us first consider (^), i.e. the conditions 



defining circular orbits. It is convenient to divide the radial potential U(u) by A(u) before differentiating with respect 
to u. Let us introduce the short-hand notation 

, ^ f 
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Eqs. (0) are then equivalent to the two equations 



aE 2 -2/3EL = -fL 2 + l, (27a) 
a' 'E 2 - 2/3' 'EL = -f'L 2 , (27b) 



where the prime denote a derivative with respect to u. Eq. ( |27b|) is homogeneous (of second degree) in E and L. We 
can therefore get a second degree equation for the ratio A = L/E. With sign conventions such that L = L z > 0, one 
finds that A is determined as 



A = A 2 (u, a) = . (28) 

7 



—2 



Inserting this solution in Eq. (27a) then determines E . Finally, we can write E and L as the following explicit 
functions of u and a: 

E = E 2 (u,a) = (a-2/3X--fX 2 )- 1/2 , (29a) 
L = L 2 (u,a) = A(a-2/3A- 7 A 2 )- 1/2 , (29b) 

where A must be replaced by Eq. (p8|). In the equations above we have added the index 2 to various functions to 
indicate that they are functions of the two variables u and a. At this stage, the results are valid for any circular orbit, 
not necessarily a member of a (specific) corotating sequence. Below, we shall determine how the spins, and thereby a, 
vary along a corotating sequence, i.e. vary as functions of u. This will give rise to functions of only one variable: u. 

From the results (|2^) one can determine the (scaled) real orbital angular momentum L, and the (scaled) real energy 
E [see Eqs. (§), (Eh and (§)]. Explicitly 



E = E 2 (u,a) = -y/ 'l + 2i/(E 2 (u,a) - 1) , (30a) 
L = L 2 (u,a) = L 2 (u,a) +a E 2 (u,a) . (30b) 

To write explicitly the condition (|2l]) let us introduce the scaled orbital angular velocity 

h = Mfl . (31) 

When a, and the spins, are fixed, Eq. (pj) gives dE Tea ] = dL so that we can compute SIq as follows from Eqs. (pC 



o / -\ dE 2 (u,a)/du 

ilo(u,a) = — — — . (32) 

d L 2 (u, a)/du 



Note that, at the LSO (corresponding to any fixed value of a), both the denominator and the numerator of Eq. Q3S 
have a simple zero, so that f^o has a well defined limit. 

The most complicated calculation is that of the spin angular velocities fii and f^- Indeed, the dependence of E rca \ 
on S a is quite involved as, besides the explicit dependence on the spins, one must also consider the implicit spin 
dependence (via Eq. (|2fj|)) of all the mass-dependent quantities: M, /i, is, a\, a 2 , . . . The simplest way to do the 
calculation is to use the general identity (24). A simplification comes from the fact that, following p2| , we can restrict 
ourselves to the symmetric case: mi = m 2 , and Si — S 2 . For such a case, it is easily seen that, when differentiating 
functions which are 1 <-> 2 symmetric (such as E Tca \) it is correct to consider quantities such as v — mim 2 /M 2 , 
<7i = 1 + 3 m2/(4mi), <j 2 = 1 + 3mi/(4m2) as fixed when varying the spins. Setting 



fti = Mfii, ai = -4 = ^a, (33) 
mi 7 



a long, but straightforward, calculation leads to the following expression for S7i : 

f>x = ^—fl^-ino^+l^/T^f^-no^ . (34) 

1 + yjl -a\ V 4 2 /16 V ^ da da J 

The functions E 2 and L 2 are that defined in Eqs. (^) above. Let us recall that Eq. ( [34]) is valid in the symmetric 
case where Si = a 2 — (8/7) a, ai = a 2 = 7/4, v = 1/4, f2i = O2, etc. . . 
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FIG. 3: Maximum binding energy configurations (LSO) obtained with various analytical and numerical methods. Empty 
(resp. filled) symbols denote irrotational (resp. corotating) systems. 'EOB 3PN' denotes the computation performed with the 
Pade approximant A(u,a 2 ) given by Eq. (H), whereas 'EOB 3PN (l-2u)' corresponds to the A (u, a 2 ) form given by Eq. (E™. 
References are as follows: Blanchet 2002 Buonanno & Damour 1999 Damour et al. 2000 [Q; Baumgarte 2000 ^8[ ; 
Cook 1994 jTfJ; Pfeiffer et al. 2000 0; Grandclement et al. Q. 



When we impose the corotating condition fii = Oo, Eq. (|34|), yields an equation to determine a\ = (8/7) a as a 
function of u. For instance, one can rewrite Eq. ( p4| ) in the form a = f(u,a). As a is a rather small quantity, we 
can solve for a by iteration: e.g., with sufficient accuracy, a ~ f(u, f(u, f(u, 0))). Having so obtained a as a function 
of u, say a — a u (u), we can then insert this result in the previous expressions to derive a "basis" of independent 
(dimensionlcss) quantities, say: 



E 
L 

mi 

TOi rr 



Ex(u) = E 2 {u,a u (u)) , 
Li{u) = L 2 (u,a u (u)) , 
ft (u,a u (u)) , 
8 - ( \ 



1+ VT 



(35) 



From this "basis" we can then compute, as functions of u, all the quantities we might need, and, in particular, the 
dimensionless quantities which are most directly derived from the numerical calculations, namely 



Mabm 


-^rcal 


M iri 




J ADM 


L + S x + S 2 


M? 

1IT 


(mf* + m 2 ") 2 




mi rr - 








mi 



1 mi 



4 mi" 

_ 1 
~ 4 



E, 



mi 
to 1 " 



[L + 2ai] 



(36) 



Characteristics of the LSO 



Figure |^ compares various analytical and numerical estimates of the maximum binding energy (.E rea i/M; rr ) — 1, 
together with the corresponding dimensionless angular velocity M- m Q, along sequences of circular orbits of equal- 
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mass black holes. Note that some sequences are corotating while others are irrotational. Several interesting lessons 
can be drawn from Fig. |[ The first one is that the analytical (EOB) results show that the effect of the corotation is 
quite small. [The reason for this smallness is discussed in detail below.] In fact, the corotation effects are smaller or 
comparable to the differences induced by making other changes in the treatment of the problem (like using a different 
resummation method for the basic EOB potential A(u), or working at the 2PN approximation instead of the 3PN 
one). 

Another important lesson drawn from Fig. || concerns the PN robustness of the EOB estimates (i.e. their stability 
under a change of PN accuracy): Indeed, we have compared the results corresponding to the 1PN approximation of 
the EOB potential A(u), namely 



A 1PN (u) = l-2u + a 2 u 2 , (37) 



A 2 ™ (u) = P\ [1 - 2u + a 2 u 2 + 2 v u 3 } , (38) 



to its 2PN approximation 



and to its 3PN approximation (jlj) (with a±(v) given by Eq. (Jl l|) ) . 

We see that, on the scale of the differences between the various estimates, and in particular, on the scale of the 
differences of the various numerical results among themselves, the EOB method yields rather close results at all 
PN approximations. This is one of the signs of the good resummation properties of the EOB method, and is to 
be contrasted with the much poorer behavior (under a change of PN accuracy) of the other analytical methods of 
estimating LSO characteristics such as: the direct use of PN-expanded equations of motion (as discussed in 
Table II of |3^]), the use of PN-expanded Hamiltonians (as discussed in |33| and |jl9|| ; see, e.g., Fig. 3 of the 
E- method (discussed in Table I, Table VII and Fig. 5 of and in Fig. 3 of the e-method (discussed in |l3| ] 

and Q) or the j-method fli) . 

For comparison purposes, we have included in Fig. ^the (irrotational and corotating) analytical predictions recently 
derived by using the the non-resummed E- method at the 3PN level Q . We note that the prediction of the non- 
resummed E-mcthod for the irrotational LSO angular velocity is 46% larger than the corresponding EOB prediction. 
The fact that the E-method yields, at 3PN, significantly larger values of (f2-Mi rr )Lso th an the EOB one, in the equal- 
mass case, is not surprising in view of the fact that it already did so in the "test mass limit" [v — ► 0), i.e. when one 
mass is much smaller than the other one. Let us introduce the dimensionless ratios, a = (OMi rr )/(f2Mi rr ) cxact , b = 
(-Ereai/Mir,. — 1) / (E lea x/Mi IT — i) cxac ( for LSO characteristics. In the test-mass limit (where the exact answers are 
known), the EOB method yields exact results, i.e. (a,b) = (1, 1) at all PN levels. By contrast, the E-method yields 
H fL4j: (8,2.914) at 1PN, (1.824, 1.315) at 2PN, and (1.275, 1.103) at 3PN. This monotone approach, from upwards, 
towards the exact answer is linked to the signs of the expansion coefficients of the function E(x). As these signs remain 
the same for all values of v one also expects the E-method to yield overestimates of f2 and E roa i in the equal-mass 
case. It is, in our opinion, an important feature of the EOB method (which is shared by none of the non-resummed 
PN approaches) that it yields exact results when v — > 0, and "robust" results (not depending very much on the PN 
accuracy used) when v is not zero, and notably when v = 1/4 (equal- mass case). On the other hand, in the equal-mass 
corotating case, the 3PN-level E-method predicts a maximal binding configuration which is within 10 % (both for f2 
and E rea i) of the 3PN EOB prediction, as well as of the HKV numerical results |l5). In view of the test-mass-limit 
behavior, and of the rather large difference in the estimates of the effects linked to corotation (see the distance between 
the diamonds in Fig. and the discussion in the next subsection), we view this 10 % agreement as accidental. 

Though the 1PN, 2PN and 3PN results are much closer in the EOB approach than in other (less resummed) PN 
approaches (compare, e.g., the results here with Fig. 3 of Ref. whose 1PN predictions would not fit on our 
Figure || above), one, however, notices that there is somewhat of a jump between the 1PN and 2PN results on one 
hand, and the 3PN ones on the other hand. This is due to the largish positive numerical value of a2 R , Eq. ([Tl]). 
Indeed, a positive 04 means a repulsive term in the radial Hamiltonian. This (relative) "repulsion" allows the LSO 
to move inwards, towards more binding. The largish positive value of a2 n (when v = 1/4) has thereby a somewhat 
significant effect on the maximal binding (towards more binding). One should, however, remember that Ref. || has 
shown that the location of the LSO was blurred by radiation damping effects. The difference between 2PN and 3PN 
predictions is probably barely observable in the GW signal emitted by coalescing black holes. [See |35| for a discussion 
of the observability of the various parameters entering the EOB framework, and for an estimate of the plausible value 
of the 4PN term af R and its effect on the LSO.l 



But the most striking feature of Fig. [3] is, on the one hand, the closeness between the EOB results and the numerical 
results of [ ^2| , and on the other hand, the huge difference between the latter results, and the numerical results of the 
conformal imaging |T(| or puncture Q IVP methods. 

To complete the information displayed in Fig. || we give in Table Q the corresponding numerical data. Table [| gives 
also data corresponding to several possible alternative versions of the basic function A(u) used in the EOB approach 
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TABLE I: Parameters of the LSO configuration according to various methods. Unless otherwise noticed, PN computations 
correspond to the EOB approach with A(u) functions given in the text and with the 3PN parameter set to 4.672, see 



Eq. ([ill). For A'" (u,a 2 ) , we have used as = 


—3. The relative 


"orbital" binding energy 


e is defined by Eq. 
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/ irr 
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0.108 


1.0019 
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2PN corot 




0.0715 


-0.0138 
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0.173 
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-0.0160 


3PN corot A(u,a 2 ) 
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-0.0168 


3PN irrot non resum. jD| 


0.129 


-0.0193 


0.786 












IVP-punct irrot Jl|] 
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(see Eqs.(15)-([T^)). In particular, note that all the Pade alternatives (A, A ,A ) yield very close results. This is 
further illustrated in Fig. [|. Even the "factorized Taylor" version (A ) (which incorporates a rather brutal way of 
forcing A(u) to have a simple zero) gives (as is clear from Fig. ||) results which are in good agreement both with other 
3PN EOB results and with the HKV ones. Though we believe that such a "factorized Taylor" version of A(u) is, a 
priori, not as good as the various Pade versions, we included it because we think that its difference with the Pade 
results gives a plausible upper limit of the "error bar" on the 3PN EOB predictions. For instance, we think that 
a plausible range for the correct maximal binding energy is E^°l° l /My TT — 1 = — 0.0157 ± 0.0011, in the corotating 
case, and E™{/M m - 1 = -0.0167 ± 0.0009 in the irrotational one. Note also that Table § (and Fig. |) include no 
data corresponding to the straightforward Taylor version ([l2|) [truncated at 0(u 5 )}. Indeed, the qualitatively different 
behavior of A Taylm (u) (without simple zero, see Fig. [j] above) prevents the presence of a minimum along the curve 
E(Q), as shown in Fig. ||. The latter Figure illustrates the fact that, even within the EOB approach, there is a contrast 
between the closeness of the Pade-type predictions and the dispersion of the Taylor-type ones (straightforward Taylor 
and factorized Taylor). Again, we can use the distance (visually displayed in Fig. ^) between Pade results and Taylor 
ones to define a (two-sided) "error bar" on the EOB Pade predictions. 



C. Spin effects in the determination of the LSO 



Before considering in more detail the characteristics of the circular orbits around the LSO, it is interesting to discuss 
the reason why there is such little difference between the irrotational and corotating cases. Indeed, if we consider, 
say, our preferred EOB potential A(u, a 2 ), at the 3PN approximation, Table | shows that the maximum corotating 
binding energy is —0.0157, while the irrotational value is —0.0167. The difference is only 1 x 10~ 3 . This difference is 
about 4 times smaller than the difference obtained in |l5| which used a non-resummed PN approach. We are going 
to see that the relatively large difference obtained in this non-resummed PN approach is rooted in the (truncated) 
perturbative nature of this approach. 

One can distinguish two leading contributions to the difference in binding energy: (i) a contribution linked to the 
(kinetic) rotational energies of the spinning black holes, and (ii) a contribution linked to the spin-orbit interaction. 
[There are also mixed terms involving, e.g., products of the spin kinetic energy with orbital effects, but one can check 
that they are numerically smaller.] In the language of the present paper, the first contribution is essentially embodied 
in the factor nii/m 1 ™ in the first Eq.([36|) above. For corotating configurations, Eq.(|34|) above approximately yields 
a-i = o-2 = 2w where u> = M m Q, ss MSI. This yields mi/m 1 " — 1 m uj 2 /2. Numerically, we have (for the irrotational 
3PN case; see Table |), u> « 0.09, so that mi/m 1 " — 1 ps 4 x 10~ 3 . This is essentially the difference "corotating - 
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irrotational" in binding energy obtained in |15[ , because, within the non-resummed PN approach of |15[ , the kinetic- 
spin contribution largely dominates over the spin-orbit interaction energy. We wish, however, to point out that there 
are consistency problems in the application of straightforward (truncated) PN-expanded methods to the determination 
of the effect of the spin-orbit interaction on the maximum binding energy. By contrast, we shall see that the EOB 
estimate of the same quantity has no consistency problems, and happens to be much larger and to nearly cancel the 
kinetic-spin contribution. 

To have a common language between the two approaches (straightforward PN approach and EOB approach), let 
us write the Hamiltonian of a binary system as H = Ho(r,L) + Hi(r, L, S), where Hi = 2LS/r 3 is the spin-orbit 
interaction to linearized order ||(|. Here, S — aiSi +(7 2 S2 denotes the effective spin (with the same notation as above 
for (7i, Si, etc.). We only consider the aligned case where the spins are parallel to the orbital angular momentum 
L. [Here, we consider the masses as given parameters, because we treat separately the effects linked to the kinetic 
energies of the spins.] To first order in Hi, the energy along sequences of circular orbits, as a function of the radial 
coordinate r, (obtained by solving dH/dr = so as to get L = L(r)), can be written as H^[r] = H^ r) [r] + H[ r) [r], 
where H^lr] is the answer corresponding to Ho(r,L) and where the first-order correction linked to Hi{r,L) ( we 
consider S as a spectator variable which will be, at the end, replaced by a function of r obtained from the corotation 
condition) reads: 



H{ r) [r] = 



„ , n d L H Q (r,L) 
Hlir ' L) ~ d? L H (r,L) drH ^ L \ 



(39) 



L=L (r) 



Here, L$(r) is the angular momentum along the radial sequences of circular orbits defined by H$(r,L). If one then 
asks (as in |[5[ M) for the variation of the energy as a function of the orbital frequency f2(r) = [d H(r, L) / 'd £]z,=.L(r)) 

one finds H^ [fi] 
correction reads: 



H^l^l] + H^'IQ], where Hq ''[*[}] is the answer corresponding to H and where the iJi-related 



r(n)r 



r (Q) 



H[ Q) [n] 



H 



(r) dH^ r) /dr 
dCl (r)/dr 



(40) 



Here, the brackets mean that one evaluates a quantity for L = L Q (r) and then changes variables via: r — r (ft). 
Having these results in hand, we can now compare the straightforward (non-resummed) PN approach to the EOB 
one. The non-resummed PN approach considers that Hi is already of high-PN order, so that, in all the 7?i-related 
corrections above one can replace the zeroth-order Hamiltonian H by its Newtonian approximation. Applying 



then the formulas above to the spin-orbit Hamiltonian (Hi cx Lr 3 ) yields the simple results: H 



(r) 



and H 



(£2) 



Hi] 



Hi], where the brackets indicate that Hi should be evaluated along the considered sequence (and 
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expressed in terms of the corresponding variable). [One easily checks that the result for H^' is equivalent to the 
formula derived in ]37j and used in Q .] 

It is interesting to note that the results are numerically different. To have such a difference is fine far from the 
LSO (and the PN calculations are both valid there), but it is problematic near the LSO. Let us first recall the 
standard approximate estimate (valid to linear order in the perturbation) of the minimum of a function of the type 
f(x) = fo(x) + fi(x) where fo(x) has a local minimum at x m { n : when working to first order in /i one easily gets 
/min ~ /o(^min) + /i(^min)- Applying this general result to the case at hand, we simply get the spin-orbit (SO) 
contributions at maximum binding by evaluating the corrections H[ r ^ or at the zeroth-order (irrotational) LSO. 
We thereby get different values (by a factor 4/3). However, it is physically (and mathematically) inconsistent to 
get different estimates of the SO-induced change in maximum binding energy. Indeed, one is considering the same 
quantity (the energy) expressed in terms of different (gauge-invariant 2 , monotonically varying) independent variables: 
r or Q. The numerical value of the minimum of this quantity should be the same (even if is reached at a different value 
of the independent variable). This consistency problem is linked to the perturbative nature of the (non-resummcd) 
PN approach which (self-consistently) treats the spin-orbit effects as being of higher PN order, and truncates away 
(most of) the "mixed" terms involving products of the PN contributions to Hq with SO effects. Note, a contrario, 
that the consistency of the non-resummed PN treatment would be recovered if one were to keep all those mixed terms. 

By contrast, the EOB approach is immune to such consistency problems. Indeed, the EOB approach does not 
approximate Hq in the correction terms by its Newtonian expression. It, instead, uses the full (non-perturbative) 
expression Eq.(p"3|). This value of Hq consistently contains an LSO, i.e. a value of r (and SI) where dH^ /dr vanishes. 
One then sees that, when evaluating (in the approximate way indicated above) the change in the maximum binding 

induced by Hi, the last correction term in h[^ vanishes, so that we get (whatever be the independent variable used) 

SEj^gQ — -ffi lso = ^1 lso - Moreover, using the full expression Eq.(|l3|) for Ho one gets the explicit expression 
(where u = M/r): 

S £ls°o B = {1 - 3 A/(2 A + ud u A)} [Hi] . (41) 

This value depends on the expression used for the EOB potential A(u), and, thereby, on the corresponding location 
of the LSO. If we consider, for simplicity, the limit v — > 0, we can use the simple form A(u) = 1 — 2u, and the 
corresponding LSO location: ulso = 1/6. This then yields SEf^o = — [-^l]- We see that the EOB approach gives a 
larger contribution than both the straightforward PN methods discussed above. This larger value is due to the fact 
that the EOB approach is non perturbative in that it treats the effective potential A(u) exactly near the LSO. 
Let us introduce the dimensionless "orbital" binding energy 

^1^ = ^-1 = 1^-1=^^-1. (42) 
MA 4 [i mi M irr v ; 

Inserting the explicit value of the spin-orbit Hamiltonian (and consistently using tlso = 6M) finally gives the following 
estimate for the spin-orbit modification of e, in the "test-mass" limit v — > 0: 

<$ SO e L so = - V3 &-*Via= -0.802 x 10~ 2 ^ 4 a . (43) 

Here, 1/4 = Av and a denotes the effective dimensionless spin parameter S/M 2 . The numerical coefficient 0.802 is 



about twice larger than the corresponding result obtained in 15 . Note that this coefficient is consistent with the 
small v limit of the result written in Eq.(4.1) of p9| |. The latter equation then shows that, in the equal-mass case 
(of interest here), v — 1/4, i.e. 1/4 = 1, the numerical coefficient giving the SO-induced change of eLSO is further 
enhanced (by a factor 1.888) to the value S so eLSO = —1.52 x 10 _2 a. Actually, the data in Table II of |2?| show that 
the latter result is valid only for values of a smaller than 0.1. For the value a ~ 0.15 relevant to us (see Table § above), 
nonlinear effects in a further enhance the change in binding energy. Interpolating between the numerical data given 
in Table II of J29| one can see that a better estimate of e for a 1=3 0.15 is eLSo ~ —0.0197. Compared to the value 
without spin-orbit contribution, eLSO = —0.0167, this yields a SO-related modification equal to S so e ~ —0.30 x 10~ 2 . 

Finally, we conclude that the EOB approach gives (at 3PN): (i) an increase of i?LSo/Mrr due to the kinetic energy 
of the spinning black holes of order +4 x 10 -3 , and a decrease linked to spin-orbit effects (as estimated within the 
EOB framework) of order —3 x 10 -3 . The net result is that -Elso /M- m increases only by +1 x 10 -3 . This analysis has 
clarified why and how the EOB approach yields a rather small final effect for the difference between the irrotational 



2 We recall that the EOB (Schwarzschild-like) radial coordinate is defined in a gauge-invariant way. 



15 




FIG. 5: Binding energy as a function of the orbital angular velocity, according to various analytical and numerical methods. 
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FIG. 6: Total angular momentum as a function of the orbital angular velocity, according to various analytical and numerical 
methods. 



case and the corotating one. We view this result as a confirmation of the generic robustness of the EOB approach. By 
contrast, we view the significantly larger change predicted by non-resummed PN approaches as further confirmations 
of the generic lack of robustness of non-resummed PN methods. An interesting lesson of this comparison is that the 
larger SO-related changes predicted by the EOB approach are crucially linked to its non perturbative nature. 



D. Evolutionary sequences 

As we said above, the LSO characteristics do not embody the really useful information contained in the various 
approximations to binary black holes dynamics. To have a complete hold on the two-body dynamics one would need 
to compare the Hamiltonians, H(r,p r , L, S a ). The EOB method provides such a complete description (which allows it 
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to describe not only the adiabatic inspiral but also the crucial transition to the plunge). However, numerical methods 
do not give access to such a multi- variable function. As a makeshift we can at least compare various functions of one 
variable linked to the two-body dynamics. In Fig. || we plot various analytical and numerical estimates of the binding 
energy of circular orbits versus the orbital frequency, while Fig. |^ plots the total angular momentum J. These figures 
show that the methods that led to good agreement for LSO characteristics in Fig. 0, also lead to good agreement for 
all circular orbits, and all available dynamical quantities. We have also included in these figures the results coming 
out of the straightforward truncation of the A-function, Eq. (|l2|). We see that it agrees rather well (compared, say, 
with the conformal- imaging or puncture data) with the helical Killing vector (HKV) data nearly up to the LSO. 

If we remember that Ref. || has shown that, a little bit above the normal (adiabatic) LSO, radiation damping 
effects dominate over the "restoring" radial potential, it is very plausible that even the pure Taylor EOB function 
A(u) would lead to a GW signal essentially indistinguishable from the one obtained from a Pade-resummed function 
A{u). Note that Figs. || andg confirm the messages of Fig. ||: the main feature is a nice agreement of all versions of 
EOB and HKV results versus a significant difference with the conformal-imaging or puncture results. When looking 
at the finer structure of the EOB predictions one sees that the 3PN predictions are much closer to the HKV data 
than the 2PN ones. 

As we said above, the multi-parameter "flexibility" of the EOB approach (a±{y), a^{y), aQ^v), . . .) can be exploited, 
to further improve, if deemed necessary, the agreement with numerical data. For instance, we have played with 
the addition of a§{v) and found that, with the Pade resummation (|l7j), a value 0,5 (|) ~ —3 further improves the 
agreement with HKV data. However, we do not consider such a best fit value of a§{y) as significant at this stage. 
Another interesting "robustness" exercise consists in modifying the general relativistic (GR) prediction ( |l0| ) for the 
3PN coefficient ( |To|) . We recall that, from the recent dimensional continuation calculation of |5|, the GR value should 
be, when v = |, et^ R = 4.672 (see Eq. (pTj)). By varying 04 (without adding any further terms a^v), ciq(v), . . . ) we 
found that a value which provides a slightly better fit is (24 (|) = 5.50. The corresponding curves are plotted in Figs. || 
and |^. On the other hand, note that the value 0,4 = 0, corresponding to the 2PN approximation only, noticeably 
worsens the agreement with the HKV data. In all cases, the reader should notice that the binding energy differences 
linked to these different choices for the 3PN contribution 04 remain small compared with the dispersion among all 
numerical results. 



IV. CONCLUSIONS 



A. Robustness of the EOB method 



The main conclusions of this work are the following. We have confirmed the robustness of the effective one-body 
(EOB) approach to binary black hole dynamics. The various ways of resumming the crucial EOB radial function 
A(u) exhibited in Figs, [l] and ^| lead to predictions for the binding energy and the angular momentum as functions 
of orbital frequency, which are very close to each other, see Figs. |, I and §, and Table J. Even the "non resummed" 
function A(u) (upper curve in Fig. [I]) leads to dynamical predictions which are quantitatively close to the "resummed" 
predictions up to the last orbits before the Last Stable (circular) Orbit (LSO) (see Figs. [3] and ^|). We have argued 
that it would probably yield indistinguishable waveforms if used, as in |6j], to study the transition between inspiral 
and plunge. [However, if one wishes, as in ||, to continue the description of the plunge nearly down to coalescence we 
expect that a non resummed A(u) will quickly lead to problems and to strong deviations from the results obtained 
from various resummed ^4(u).] Another robustness feature of the EOB approach is its good post-Newtonian (PN) 
convergence properties. This was illustrated in Figs. ||, || and || where the 1PN, 2PN and 3PN predictions are seen 
to be all much closer to the helical Killing vector (HKV) data than to the other numerical data. Note, however, 
that the best agreement with HKV data is obtained for the 3PN prediction. We emphasized also the consistency 
(and robustness) of the EOB approach in the treatment of the changes induced by spin effects. This contrasts with 
non-resummed PN approaches which predict different changes when one uses different parameters to run along a 
corotating sequence. 



B. Extension of the EOB method to corotating systems 



The most novel aspect of the present work is that we have shown how to apply to corotating sequences the recently 
defined |2^] extension of the EOB approach to spinning black holes. The effects linked to the spins turn out, finally, 
to be rather small. For instance, they are significantly smaller than the effect of using the 3PN-accurate EOB 
Hamiltonian instead of a 2PN-accurate one, see Fig. 0. We discussed the reasons behind this behavior (essentially a 
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near cancellation between an energy increase linked to the spin kinetic energy and an energy decrease linked to the 
rather strong spin-orbit interaction predicted by the EOB method). 

C. Good agreement between EOB and HKV results 

Our findings fully confirm what was announced in [^2| : the new HKV numerical data are much closer to analytical 
results than previous numerical data. At face value, this is a very encouraging fact which suggests that one can 
now implement the philosophy which was at the basis of the EOB approach ||: to use a (resummed) analytical 
approach to describe binary black holes not only during early inspiral (which was previously thought to be the only 
possible domain of validity for an analytical approach), but also during late inspiral and, most importantly, during 
the transition between inspiral and plunge. Numerical relativity techniques can then be used only for describing the 
rest of the plunge and the coalescence. [As noted in [^0) one might also simplify matters by describing the end of the 
coalescence by the "close limit approximation" .] A crucial requirement for successfully merging together analytical 
and numerical results is to be able to match the analytical description just after the last stable orbit (LSO) with 
suitable numerical initial data. This paper is an encouraging step in this direction, in view of the robust and close 
agreement between EOB results and HKV ones. 

However, this agreement raises many questions. A first question is linked to the "approximations" used in the present 
implementation of the HKV philosophy. Indeed, though this method is not, in principle, limited to conformally flat 
data, its present implementation has, for simplicity, imposed conformal flatness of the spatial metric, and has reduced 
the ten equations (for ten unknowns) to solve to a subset of five equations (for five unknowns). The reason why 
this restriction is problematic is that, from the analytical side, the values of the successive coefficients in the various 
functions entering the EOB formalism depend quite sensitively on the exact form of the metric, and, in particular, 
on the fact that it is not conformally flat. For instance, it was shown in jl4|] that both the "effective potential" 
A(u), and the energy map ([!]), are modified, at the 2PN approximation, if one truncates Einstein equations to impose 
conformal flatness. The corresponding modification at the 3PN approximation is currently not known. However, 
these modifications might turn out not to be critical in view of the robustness of the EOB predictions. Fig. || shows 
that even if we use the EOB Hamiltonian at the 1PN approximation (which, in fact, essentially coincides with the 
"Newtonian" approximation, i.e. the EOB reformulation of the "test- mass approximation"), one is much closer to 
HKV data than to the other numerical data. 



D. Plausible explanation of the disagreement with IVP results 

The most important question posed by the agreement analytical/HKV is to understand why the IVP approaches 
give (consistently among themselves) drastically different results. After all, both types of numerical approaches use 
(at present) the same simplification of spatial conformal flatness (7^ = ip 4 Sij), and they both numerically construct 
solutions to the Hamiltonian and momentum constraints (approximate solution to the momentum constraint for 
HKV). Seen from this technical point of view, the difference between the two schemes lies in the way of solving the 
momentum constraint. In the IVP method |l^, [l^, ^8), it is solved by means of a simple ansatz (Bowen-York extrinsic 
curvature f3§| or its isometric version fl39| ) for the second fundamental form whose weak- field limit is reasonable, 
while in the HKV method [^l], |2^] one determines an essentially unique solution for Kij by the requirement that it 
be induced by restricting a helically-symmetric spacetime to a Cauchy hypersurface (see also |2j| for a discussion). 
We think that the basic physical difference between the two approaches, which underlies their drastically different 
predictions, is the following. Though it must be admitted that the imposition of a helical symmetry can only be 
an approximation, it is a physically well motivated approximation which has the great virtue of selecting a specific 
solution of the constraints. This is very similar to what happens in the post-Newtonian schemes, especially in their 
ADM version. Indeed, in the Hamiltonian ADM approach one uses (when dealing with the near-zone field) the 
post-Newtonian ansatz (c -1 dt h <C d x h) to derive and solve an iterative set of coupled elliptic equations for 7y 
and K^. This selects a specific solution of the constraints. Intuitively, one can think of this (essentially) unique 
solution selected by the PN method as the only solution of the constraints which contains the "right" amount of free 
gravitational wave (GW) data (hf^, 7r^ T ) which has been generated by the many previous orbits of the system, and 
which has therefore reached a quasi-stationary state. Any other amount of GW data would be a priori allowed at 
a specific moment of time but would not correspond to a quasi-stationary state, and would therefore be expected 
to "fly away" with the velocity of light, i.e. to evolve in a violently non quasi-stationary way, with c^ 1 dfh ~ d x h 
instead of c _1 dth <C d x h. We think that this PN understanding of the selection, by reduction to elliptic equations, 
of a specific "quasi-stationary" solution applies, mutatis mutandis, to the HKV approach. The imposition of a helical 
Killing vector selects the only "quasi-stationary" solution corresponding to a steady situation generated by a slow 
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inspiral, while the arbitrary choice, in the other methods, of a technically simple solution of the momentum constraint 
corresponds to a drastically non-steady state, containing a wrong amount of free GW data 3 , which is expected to "fly 
away" as soon as one tries to evolve the system in time. If this view is correct, it diminishes the signification of the 
"effective potential" approach applied to the initial data of |n| [l?], [p| . Indeed, these data correspond to a non-steady 
state to which it is incorrect to apply any energy minimization principle. 



E. Fitting the EOB Hamiltonian to numerical data 

We have started to exploit the "flexibility" of the EOB approach (already emphasized in ^9|), i.e. the use of the 
expansion coefficients appearing in the EOB approach as parameters to be fitted. For instance, Figs. |] and ^| exhibit 
the dynamical predictions obtained from a modified 3PN coefficient, <24 estfit (j) ~ 5.50 instead of R (j) ~ 4.67. 
The modified value of leads to a better agreement with the HKV data. As we are not adding any higher PN 
contributions, a^,(v)u b + ag(^)u 6 + . . . , such a best fit value for can be viewed as a way of mimicking the effect of 
such (missing) higher PN terms. 

We resisted the temptation to (introduce and) vary more a n {v) parameters up to getting a really close agreement 
with HKV data, because, at this stage, such a formal exercise would not be justified. It is, however, important to 
keep in mind the suggestion made in [f29| that a suitable "numerically fitted" EOB Hamiltonian may be a very useful 
tool for combining the advantages of analytical and numerical methods. If we look ahead to the problem of exploring 
the dynamics of two arbitrarily spinning black holes, one will have to face an extremely large parameter space. It 
seems clear that numerical relativity will not be able (before many years) to cover densely this parameter space. On 
the other hand, one can use some sparse numerical data to determine free parameters introduced in a generalized 
EOB Hamiltonian. Then this "numerically fitted" EOB Hamiltonian can be used to interpolate between the sparse 
numerical data. 
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